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We  have  used  large-eddy  simulation  with  an  immersed  boundary  method  to  study 
turbulent  flows  over  distributions  of  uniform  height,  staggered  cubes.  The  computational 
domains  were  designed  such  that  both  the  roughness  sublayer  and  a  region  of  the  inertial 
layer  are  resolved.  With  this,  we  record  vertical  profiles  of  time  series  of  fluctuating 
streamwise  and  vertical  velocity  at  different  locations  throughout  the  domain.  Contour 
images  of  these  fluctuating  quantities  shown  relative  to  elevation  and  time  are  studied; 
contour  images  of  Reynolds  shear  stresses  owing  to  ‘sweeps’  and  ‘ejections’  are  also 
studied.  These  images  show  that  periods  of  momentum  excess  (deficit)  in  the  inertial- 
layer  precede  excitation  (subdual)  of  cube-scale  coherent  vortices  in  the  roughness 
sublayer.  We  compute  this  time  lag  (termed  advective  lag)  and  demonstrate  that  it  scales 
linearly  with  wall-normal  elevation.  The  advective  lag  is  attributed  to  coherent,  low- 
and  high-momentum  regions  in  the  aloft  inertial  layer.  Vortex  identification  is  used  to 
illustrate  the  presence  of  hairpin  packets  encapsulating  low-momentum  regions.  Based 
on  this,  the  reported  inclination  angle  associated  with  hairpin  packets  is  used  to  guide  the 
development  of  a  model  for  prediction  of  advective  lag  with  height.  The  model  captures 
the  advective  lag  profiles  reasonably  well.  In  the  interest  of  generality,  additional  cases 
of  flow  over  homogeneous  roughness  (aerodynamic  drag  imposed  with  the  equilibrium 
logarithmic  law)  are  considered.  We  again  observe  that  advective  lag  scales  linearly 
with  wall-normal  elevation.  Advective  lag  predictions  from  the  aforementioned  model 
agree  well  with  results  for  these  cases. 

Keywords:  large-eddy  simulation;  roughness  sublayer;  atmospheric  surface  layer; 
coherent  motions;  hairpin  packets;  advective  lag 


1.  Introduction  and  motivation 

The  dynamics  of  turbulent  flows  over  rough  walls  are  important  in  numerous  settings. 
The  performance  of  vapour  power  systems,[l]  the  efficiency  of  naval  vessels, [2]  and  the 
economics  of  pipeline  transport  are  the  examples  of  contemporary  applications  greatly 
influenced  by  turbulent  momentum  exchanges  in  close  proximity  to  the  wall.  Meanwhile, 
turbulent  momentum  transport  in  atmospheric  boundary  layer  flows  is  of  critical  importance 
to  the  performance  of  modern  wind  farms, [3]  aerodynamics  of  vegetative  canopies  [4,5]  and 
urban  environments, [6-9]  and  geomorphological  processes  associated  with  the  evolution  of 
aeolian  desert  landscapes. [10]  The  aforementioned  engineering  and  geophysical  examples 
are  distinctly  different  to  smooth  wall  turbulence,  owing  to  the  presence  of  a  distribution 
of  obstacles  of  height,  h,  that  protrude  into  the  inertial  layer  of  the  flow.  These  obstacles 
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absorb  momentum  through  pressure  drag,  induce  flow  separation,  and  serve  to  produce 
obstacle-scale  coherent  motions  that  occupy  the  region  between  the  wall  and  approximately 
three  obstacle  heights  above  the  roughness  sublayer.  [1 1,12]  In  the  roughness  sublayer,  the 
aerodynamic  drag  distribution  due  to  the  presence  of  obstacles  results  in  an  inflected  mean 
streamwise  velocity  profile  at  the  approximate  average  obstacle  height. [13, 14]  As  a  result, 
the  mean  flow  gradient  exhibits  its  maxima  at  the  inflection  [10]  (not  the  wall),  and  the  tur¬ 
bulence  kinematics  are  fundamentally  different.  For  flow  over  vegetative  canopies,  Raupach 
et  al.  [15]  demonstrated  that  flow  in  the  roughness  sublayer  resembled  a  turbulent  mixing 
layer  with  positively  and  negatively  skewed  streamwise  and  vertical  velocity  fluctuations, 
respectively,  and  Reynolds  stresses  composed  predominately  of  ‘sweeps’.  With  this,  the 
turbulence  structure  is  characterised  by  Kelvin-Helmholtz  spanwise  ‘rollers’,  which  origi¬ 
nate  at  the  velocity  profile  inflection  and  undergo  a  downstream  transformation  leading  to 
hairpin  packets  [4]  (while  remaining  contained  within  the  roughness  sublayer).  Recently, 
Ghisalberti  [16]  demonstrated  the  existence  of  an  universality  in  roughness  sublayer  statis¬ 
tics  for  flows  over  diverse  canopies,  and  introduced  the  term  ‘  obstructed  shear  flow’  to 
categorise  such  flows.  Within  the  sublayer,  the  turbulence  geometric  macroscale  is  set  by 
the  vorticity  thickness, [15, 17]  <$w  =  U(h)/(dU(z)/dz)\h,  and  streamwise  spacing  of  vortex 
cores,  Ax  =  2n  L^,  where  U(z)  and  L}jj  is  the  mean  streamwise  velocity  and  the  integral 
length  scale,  respectively.1  The  mixing  efficiency  associated  with  these  Kelvin-Helmholtz 
roughness  sublayer  motions  -  which  is  commonly  quantified  via  the  correlation  coefficient, 
ruw  =  u2t /MrmsWms,  where  ur  is  friction  velocity  and  rms  denotes  root  mean  square  -  ex¬ 
ceeds  the  value  exhibited  by  a  logarithmic  boundary  layer  by  approximately  50%.  [16, 18] 
Thus,  turbulence  in  the  roughness  sublayer  is  characterised  by  vigorous  mixing  and  complex 
structural  attributes. 

Above  the  roughness  sublayer  -  in  the  inertial  layer  -  the  mean  flow  profile  exhibits 
logarithmic  scaling  and  the  turbulence  structure  resembles  smooth  wall  flow.[8,13,14]  That 
is,  the  domain  is  occupied  by  persistent  streamwise-elongated  coherent  parcels  of  relatively 
low  and  high  momentum.  For  smooth  walls,  such  low-momentum  regions  (LMRs)  mean¬ 
der^  19]  and  are  encapsulated  by  hairpin  packets  [20-27]  at  the  interface  between  zones 
of  quasi-uniform  momentum.  [28]  For  rough  walls,  the  presence  of  persistent  structures 
is  also  reported,[8,25]  however,  experiments  have  indicated  that  A-scale  coherent  motions 
associated  with  the  roughness  sublayer  seemingly  reduce  the  lengths  of  logarithmic-layer 
coherent  motions. [25,29]  Recently,  Coceal  et  al.  [8]  used  direct  numerical  simulation  (DNS) 
to  study  flow  over  a  staggered  array  of  uniform  height  cubes  [6]  with  characteristic  scale, 
h/H  =  i,  where  H  is  the  boundary  layer  depth.  They  illustrated  the  existence  of  hairpin 
packets  (and  ‘cane’  structures,  which  are  inclined  coherent  parcels  with  only  one  leg  of 
the  hairpin  [24]  around  the  LMR).  Moreover,  they  used  approximated  conditional  averages 
[20,30]  to  illustrate  the  significant  reductions  in  streamwise  coherence  of  their  rough  wall 
turbulence  relative  to  a  smooth  wall.  These  efforts  to  characterise  the  structural  attributes 
of  smooth  and  rough  wall  turbulence  have  typically  focused  on  spatial  characteristics. 
Furthermore,  most  previous  studies  address  either  instantaneous,  ensemble-averaged  or 
time-averaged  statistically  stationary  turbulence  statistics.  Here,  we  have  used  large-eddy 
simulation  (LES)  to  study  the  temporal  dynamics  of  turbulent  flows  over  such  roughness. 
This  is  accomplished  by  recording  time-series  flow  statistics  across  the  depth  of  the  bound¬ 
ary  layer  at  a  variety  of  horizontal  positions  (thus,  we  use  LES  to  generate  data-sets  that 
would  more  typically  be  associated  with  field  campaign  data  from  a  tower  equipped  with 
sonic  anemometers).  [3 1]  An  advantage  of  our  approach  is  that  it  would  provide  insight  into 
the  interpretation  of  laboratory  or  field  experimental  results,  which  are  often  collected  as 
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Figure  1.  Illustration  of  [6]  block  case  (‘C20S’)  considered  as  lower  topography  in  this  study:  (a) 
perspective  image  showing  freestream  flow  direction,  Uo  and  (b)  planform  image  with  indication  of 
positions  at  which  time-height  velocity  data  are  recorded.  Owing  to  the  topography  attributes,  the 
four  points  effectively  capture  temporal  dynamics  at  all  points  in  the  domain.  The  above  topography 
is  labelled  here  Case  SC2. 


time  series.  Moreover,  illustrations  of  ‘time-height’  colour  floods  of  fluctuating  stream- 
wise  velocity  offer  a  novel  perspective  on  the  presence  of  the  aforementioned  coherent  flow 
structures  and  clearly  show  the  extent  to  which  activity  in  the  roughness  sublayer  is  influ¬ 
enced  by  inertial-layer  dynamics.  Moreover,  by  incorporating  these  results  with  arguments 
relating  to  structural  attributes  of  coherent  flow  structures  in  wall-bounded  flows,  we  have 
been  able  to  develop  a  predictive  model  for  the  reported  advective  lag.  We  stress  that  the 
results  presented  here  (and  related  predictive  capability)  are  not  expected  to  apply  only  to 
flow  over  urban-like  complex  topographies.  Thus,  the  predictive  capability  could  be  used, 
for  example,  in  conjunction  with  an  atmospheric  LiDAR  system  for  control  of  wind  farms. 


1.1.  Present  study 

We  use  LES  to  model  flow  over  topographies  composed  of  a  staggered  distribution  of 
wall-mounted  cubes  (one  of  the  cases  considered  is  ‘C20S’  from  the  wind  tunnel  study  by 
Cheng  and  Castro  [6],  and  considered  in  the  more  recent  DNS  study  by  Coceal  et  al.  [8]). 
Figure  1(a)  shows  a  perspective  view  of  the  array  of  cubes  with  edge  length,  h.  Figure  1(b) 
shows  the  arrangement  in  planform,  with  indication  of  the  streamwise  spacing,  s/h ,  between 
rows  of  cubes  (s/h  =  1  for  all  cases  in  this  study).  Figure  1(b)  includes  points  xi,x2,  X3,  and 
X4 ;  during  LES,  time  series  of  streamwise  and  vertical  velocity  have  been  recorded  across 
the  depth  of  the  boundary  layer  at  these  points.  This  is  better  illustrated  in  Figure  2(a), 
which  shows  the  cubes  and  streamwise  spacing  (for  the  cube  cases,  they  coordinates  of 
positions  X1-X4  are  equal  and  set  to  intersect  the  centre  of  cubes;  the  x  coordinates  are 
varied  such  that  xi  is  precisely  at  the  centre  of  a  cube,  x2  is  centred  between  the  upwind 
and  downwind  row  of  cubes,  X3  is  precisely  at  the  centre  of  adjacent  cubes,  and  X4  is 
centred  between  the  upwind  and  downwind  row  of  cubes).  For  generality,  we  also  model 
flow  over  homogeneous  roughness  lengths,  shown  in  Figure  2(b).  The  novelty  of  this  work 
resides  in  the  insights  gained  from  considering  time-height  contours  of  fluctuating  velocity 
components  in  the  context  of  the  aforementioned  (and  comprehensively  studied)  structural 
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(a)  x\  x2  x3  X4 


(b)  xc 


Figure  2.  Illustration  of  computational  domains  considered  in  this  study:  (a)  uniform  height,  stag¬ 
gered  array  of  blocks  [6]  with  h  varying,  s/h  =  1  and  (b)  homogeneous  roughness,  r0.  For  panel 
(a),  see  also  Figure  1.  For  all  cases,  H  =  1000  m  and  other  simulation  parameters  are  summarised 
in  Table  1 .  Illustrative  locations  of  positions  at  which  time  series  of  streamwise  and  vertical  veloc¬ 
ity  are  recorded  are  shown  for  discussion  (vertical  red  lines).  In  both,  an  illustrative  time-averaged 
streamwise  velocity  profile  is  shown,  where  (7o  is  the  ‘outer’  freestream  velocity. 


attributes  [8,20-22]  of  turbulence  in  the  roughness  sublayer  [8,25]  and  inertial  layer.  The 
results  (Sections  3. 1-3.3)  and  scaling  arguments  (Section  3.4)  demonstrate  the  influence 
of  inertial-layer  structures  on  dynamics  of  the  roughness  sublayer.  Section  2  contains  a 
description  of  the  LES  code  and  validation  against  an  available  wind  tunnel  data-set  [6];  a 
description  of  the  cases  considered  for  this  study  is  also  provided  in  Sections  2.  Section  3 
shows  a  series  of  results  that  lead  to  a  semi-empirical  model  for  reported  temporal  dynamics 
of  the  urban  sublayer.  The  influence  of  resolution  (spatial  and  temporal)  is  tested  in  Section 
3.5;  we  report  that  advective  lag  does  not  exhibit  dependence  on  resolution.  Concluding 
remarks  are  provided  in  Section  4. 


2.  Method  and  cases 

We  consider  flow  over  Figure  2(a)  and  2(b)  topographies  with  attributes  summarised  in 
Table  1.  During  LES,  the  spatially  filtered  three-dimensional  incompressible  momentum 
transport  equations  are  solved  at  high-Reynolds  number  [32-35]  for  a  neutrally  stratified 
(i.e.  no  buoyancy  forces)  turbulent  boundary  layer  without  Coriolis  accelerations: 

du  1  1  1 

- 1 — V  (u  ■  u)  —  u  x  co  = - V p  —  V-r  +  II  H — /,  (1) 

dt  2  p  p 

where  p  is  the  modified  pressure,  r  is  the  subgrid-scale  (SGS)  stress  tensor,  and  n  = 
{u2JH,  0,  0}  is  the  mean  pressure  gradient  in  the  streamwise  direction,  where  uT  is  the  shear 
velocity  (ur  =  0.45  m/s)  and  H  is  the  boundary  layer  depth  (II  —  1000  m).  A  solenoidal 


Table  1.  Summary  of  simulation  parameters  and  domain  attributes.  For  all  cases,  H  =  1000  m. 


Description 

Name 

L/H 

h/H 

zJH 

Nx  x  Nv  x  Nz 

Staggered  cubes 

SCI 

4 

1/8 

2  x  10-5 

128  x  128  x  128 

Staggered  cubes 

SC2 

2 

1/4 

2  x  10~5 

128  x  128  x  128 

Uniform  roughness 

U1 

4 

- 

1  x  10-3 

128  x  128  x  128 

Uniform  roughness 

U2 

4 

- 

1  x  10"2 

64  x  64  x  64 

Uniform  roughness 

U3 

4 

- 

1  x  10“3 

64  x  64  x  64 

Uniform  roughness 

U4 

4 

- 

1  x  10"4 

64  x  64  x  64 

Uniform  roughness 

U5 

4 

- 

1  x  10”5 

64  x  64  x  64 

Uniform  roughness 

U6 

4 

- 

1  x  10~6 

64  x  64  x  64 
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velocity  field  is  maintained  by  computing  the  divergence  of  Equation  (1),  applying  the 
incompressibility  condition,  V  ■  u  =  0,  and  solving  the  resulting  pressure  Poisson  equation 
for  a  pressure  correction.  Note  also  that  the  viscous  stresses,  vV2u,  are  omitted  from 
Equation  (1),  owing  to  the  high-Reynolds  number  typical  of  atmospheric  surface  layer 
flows  (the  ‘macroscale’  Reynolds  number  is  Rer  =  UqH/v  ~  (?(109)).  The  deviatoric 
component  of  r,y  is  evaluated  using  the  eddy  viscosity  modelling  approach, 


r  —  -<5Tr(r)  =  —2  v,S, 


(2) 


where  v,  =  (Cs  A)2|S|  is  the  turbulent  viscosity,  Cs  is  the  Smagorinsky  coefficient, [36]  A 
is  the  filter  size,  S  =  j(Vm  +  ViiTr)  is  the  resolved  strain-rate  tensor,  and  |5|  =  (2 S:S)1/2 
is  magnitude  of  the  resolved  strain-rate  tensor.  For  this  work,  Cs  is  evaluated  dynam¬ 
ically  during  LES  with  the  Lagrangian  scale-dependent  dynamic  SGS  model  of  Bou- 
Zeid  et  al.[34].  Pseudospectral  discretisation  is  used  in  the  horizontal  directions,  while 
vertical  gradients  are  evaluated  with  centred  second-order  finite  differencing.  Periodic 
boundary  conditions  are  imposed  on  the  vertical  planes  of  the  domain,  owing  to  spec¬ 
tral  discretisation  in  the  horizontal  directions.  At  the  domain  top,  the  zero-stress  Neu¬ 
mann  boundary  condition  is  imposed  on  streamwise  and  spanwise  velocity,  du/dz\z/H=o  = 
dv/dz\z/H=o  —  0.  The  zero  vertical  velocity  condition  is  imposed  at  the  domain  top  and 
bottom,  w(x,  y,  z/H  =  0)  =  w(z,  y,  z/H  =  1)  =  0.  Zero-stress  Neumann  boundary  con¬ 
ditions  are  imposed  on  the  pressure  Poisson  equation  solution  at  the  domain  top  and  bot¬ 
tom,  dp/dz\z/H=o  =  dp/dz\z/H=i  =  0.  The  Adams-Bashforth  time-advancement  scheme 
is  used  for  temporal  integration  of  Equation  (1).  The  nonlinear  advection  term  is  de-aliased 
in  Fourier  space  with  the  3/2  rule  [37];  this  is  necessary  since  aliasing  errors  may  contami¬ 
nate  the  smallest  resolved  scales  of  the  flow,  compromising  predictions  of  the  SGS  models. 
The  computational  domain  discretisation  is  A.v  =  Ay  =  Lx/Nx ,  and  Az  =  H/Nz,  where  L 
is  indicated  in  Table  1  and  Nx  —  Ny  =  Nz.  The  computational  domain  is  staggered  in  the 
vertical  direction;  the  first  computational  level  for  u  and  v  is  located  at  elevation  \AZ. 

Flere,  we  consider  cubic  roughness  cases  with  h/H  —  1/4  and  1/8  (SCI  and  SC2).  For 
cases  SCI  and  SC2,  the  cube  spacing  is  a  block  height  (i.e.  s/h  =  1).  We  also  consider 
six  cases  of  homogeneous  momentum  roughness,  zo  (U1  to  U6).  The  computational  mesh 
for  all  cube  simulations  is  discretised  with  Nx  =  Nv  =  Nz  =  128,  while  the  homogeneous 
roughness  cases  are  adequately  resolved  with  Nx  =  Ny  =  Nz  =  64.  [34]  Table  1  summarises 
the  simulation  attributes,  where  L/H  is  the  domain  horizontal  length  (streamwise  and 
spanwise  extent  equal  for  all  cases).  We  record  time  series  of  streamwise,  u/uT,  and  vertical, 
w/uT,  velocity  at  all  vertical  computational  levels  across  the  depth  of  the  boundary  layer  at 
positions  x\  to  X4  for  the  cube  cases  (Figure  1(b)).  Since  the  urban  topography  cube  height 
is  uniform,  selecting  positions  x\  to  X4  ensures  that  we  capture  representative  data  at  all 
other  ‘cell-centred’  positions.  For  the  homogeneous  roughness  cases,  we  record  time-series 
data  at  the  geometric  centre  of  the  domain. 

The  immersed  boundary  method  (IBM)  used  here  is  based  on  a  wall  model,  which 
depletes  momentum  based  on  the  unit-normal  area  on  which  the  flow  impinges  [33,38]  and 
the  associated  drag  is  added  as  a  body  force  in  the  momentum  transport  equation.  For  cases 
U 1  to  U6  (homogeneous  rough),  the  equilibrium  logarithmic  law  [39,40]  is  used  exclusively 
to  impose  aerodynamic  drag: 


TT  12  C 

kU  ut 
In  (z/z0) .  U’ 


(3) 
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where  i  =  1  and  2  corresponds  with  the  streamwise  and  spanwise,  respectively,  k  is  the  von 
Karman  constant  (k  —  0.4),  zq  is  a  (momentum)  roughness  length,  m,  and  w,  denote  grid-  and 
test-filtered  quantities,  respectively,  and  U(x ,  y)  =  ,  y)  +  v  (x,  y)j 1  ~  is  magnitude  of 

the  local  test-filtered  velocity.  See  Table  1  for  the  particular  roughness  lengths  used.  Here,  we 
follow  Bou-Zeid  et  al.  [34]  by  using  the  test-filtered  velocities  for  computing  the  wall  stress 
(Equation  (3)).  This  approach  is  typically  used  for  modelling  flows  over  heterogeneous 
[41]  or  complex  [33]  topographies,  since  it  serves  to  reduce  variance  of  the  streamwise 
and  spanwise  velocity  components  close  to  the  wall,  thereby  improving  prediction  of  the 
logarithmic  law.  For  cases  SCI  and  SC2  (Table  1),  we  use  the  equilibrium  logarithmic  law 
for  regions  with  h(x,  y)  =  0  (Equation  (3))  in  conjunction  with  an  IBM  for  regions  with 
h(x,  y)  >  0.  The  present  IBM  technique  imposes  aerodynamic  drag  owing  to  the  obstacles 
by  computing  a  body  force,  /,  in  Equation  (1)  due  to  spatial  variation  of  the  topography, 
h(x,  y).  For  regions  where  h(x,  y)  >  AJ2,  the  orographic  drag  force  imposed  on  the  flow 
by  h(x,  y)  is 


F  =  -  j  pwndS,  (4) 

where  «,  is  the  unit  normal  vector  to  h(x,  y)  and  pw  is  the  resolved  wall  pressure  acting  on 
h(x,  y).  This  IBM  technique  (previously  named  the  gradient-based  drag  modelling  tech¬ 
nique,  Anderson  and  Meneveau  [33],  and  later,  Anderson  [38])  is  used  to  model  pressure 
drag  forces  due  to  pw.  After  division  by  density,  p,  and  local  computational  cell  volume, 
AxAyAz,  Equation  (4)  reduces  to  a  drag  force  per  unit  mass  required  in  Equation  (1): 

-  /  = - i-  f  PwndS  *  -&R(fi  ■  V/i)  J-,  (5) 

P  P A2AZ  Js  A, 

where  R(x)  is  the  ramp  function  (R(x)  =  x  if  x  >  0,  and  R(x)  —  0  if  x  <  0)  and  dh/dxk  ( k  = 
1  and  2)  is  the  gradient  of  h(x,  y)  in  the  x  and  y  directions.  This  use  of  the  ramp  function 
isolates  frontal  areas  of  h(x,  y)  on  which  u,  impinges.  We  assume  the  drag  coefficient,  Q 
=  2,  for  all  cases,  and  therefore  the  typical  \  C,i  factor  is  omitted  in  Equation  (5)  (this 
implies  complete  depletion  of  incoming  momentum).  The  approach  was  tested  against 
numerous  data-sets  available  in  the  literature  for  flow  over  different  kinds  of  topography  - 
blocks,  sinusoids,  ellipsoidal  mounds  -  and  in  all  cases  the  performance  was  satisfactory 
(agreement  of  time-  and  plane-averaged  streamwise  velocity  profiles  within  10%).  The 
Anderson  and  Meneveau  [33]  technique  is  somewhat  unique,  as  it  is  used  to  specify  drag 
imposed  by  topography  resolved  in  the  horizontal  directions  but  not  the  vertical  (i.e.  h(x, 
y)  does  not  exceed  the  height  of  the  first  computational  mesh  level,  A./2).  Here  we  use 
the  vertically  resolving  version  of  this  modelling  approach, [38]  where  Equation  (5)  models 
drag  associated  with  ‘cut’  cells  while  the  velocity  at  fully  immersed,  internal  cells  is  set 
to  zero.  For  comparison  of  time-  and  plane-averaged  flow  statistics  from  this  modelling 
technique  against  literature  data-sets  for  a  variety  of  topographies,  the  interested  reader 
may  consult  Anderson  [38]. 

All  velocities  are  normalised  by  shear  velocity,  uT.  For  cases  U1  to  U6,  shear  velocity  is 
determined  simply  as  uT  —  (r^/p)  I/“;  for  cases  SCI  and  SC2,  shear  velocity  is  determined 
based  on  the  maximum  Reynolds  shearing  stress  occurring  in  the  roughness  sublayer: 

uT  =  max([\{u'w')Xiyj  +  {xxz)  x,y_t\\12), 


(6) 
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Figure  3.  Vertical  profiles  of  time-  and  plane-averaged  streamwise  velocity  (solid:  present  LES; 
dashed:  direct  numerical  simulation  by  Coceal  et  al.  [8])  in  linear-linear  axis  scaling  (a)  and 
logarithmic-linear  axis  scaling  (b).  In  panel  (b),  a  logarithmic  profile  is  included  (see  figure  an¬ 
notation)  for  z0/H  =  1CV2.  Vertical  profiles  of  time-averaged  streamwise  velocity:  (c)  position  xt 
in  Figure  1(b);  (d)  position  x2  in  Figure  1(b);  (e)  position  x3  in  Figure  1(b);  and  (f)  position  x4  in 
Figure  1(b).  In  panels  (c)  to  (f),  solid  line  and  circles  denote  results  from  present  LES  and  experiments 
(Cheng  and  Castro  [6],  case  ‘C20S’),  respectively. 


where  shearing  stress  associated  with  resolved  fluctuations  is  evaluated  under  the  horizontal 
statistical  heterogeneity  condition, [3, 5]  u'  —  u  —  {u)x,y,i,  and  (. .  .)a  denotes  averaging 
over  dimension  a.  In  spite  of  concerns  posed  by  this  approach  owing  to  mean  flow  spatial 
gradients  in  the  sublayer,  we  have  used  this  approach  and  accomplished  close  agreement 
with  the  literature  data-sets,  as  shown  in  Figure  3.  We  stress  that  for  cases  SCI  and  SC2, 
zq  is  based  on  regions  with  h(x,  y)  —  0;  additional  aerodynamic  drag  due  to  the  obstacles 
( h(x ,  y )  >  0)  is  imposed  with  the  IBM.  In  Section  3.2,  discussion  regarding  the  ‘effective’ 
roughness  length,  no,  Eff.,  for  cases  SC  1  and  SC2  is  provided  (also  summarised  in  Table  2). 

In  order  to  demonstrate  fidelity  of  the  IBM  [38]  technique  used  within  the  present 
LES  (Equation  (1))  for  resolving  sublayer  turbulence,  results  of  a  validation  case  are 
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Table  2.  Summary  of  simulation  results,  zo,  e ttJH  is  evalu¬ 
ated  by  fitting  a  logarithmic  velocity  profile  with  varying  z0 
to  plane-averaged  streamwise  velocity  profiles  from  LES. 


Name 

Zo,  E sJH 

SCI 

5.412  x  10"3 

0.130 

SC2 

1.08  x  10-2 

0.254 

U1 

m 

1 

O 

X 

0.012 

U2 

1  x  10-2 

0.023 

U3 

m 

1 

O 

X 

0.023 

U4 

1 

0 

X 

0.023 

U5 

1  X  10“5 

0.023 

U6 

X 

0 

1 

OS 

0.023 

presented  here.  The  case  is  in  fact  SC2  from  Table  1;  the  literature  data-sets  are  sourced 
from  Cheng  and  Castro  [6]  (wind  tunnel)  and  Coceal  et  al.  [8]  (DNS).  Although  the 
present  LES  Reynolds  number  is  much  larger  than  for  the  experimental  and  DNS  data¬ 
sets,  all  the  below  comparison  is  performed  with  outer-scaled  statistics.  Since  outer-scaled 
turbulence  statistics  exhibit  Reynolds  number  dynamic  similarity  (‘fully  rough’  conditions 
[42]),  meaningful  comparison  between  the  present  LES  results  and  literature  data-sets 
from  DNS  or  experiments  can  be  made  (Figure  3).  Figure  1(b)  shows  a  plan  view  of  the 
topography,  where  the  flow  direction  is  aligned  in  the  x-direction.  Figure  1  shows  locations 
xi,x2,x3,  and  X4 ,  at  which  vertical  profiles  of  time-averaged  streamwise  velocity,  ( u),/uz , 
are  provided  in  the  literature  data-sets. [6,8]  Figure  3  shows  vertical  profiles  of  {u),/uz  from 
the  present  LES  (solid  line),  along  with  experimental  [6]  (circles)  and  DNS  [8]  (dashed 
line).  Figure  3(b)  shows  Figure  3(a)  profile,  plotted  in  logarithmic-linear  axis  scaling;  for 
discussion,  a  logarithmic  profile  has  been  added,  to  illustrate  that  the  LES  profile  (solid, 
curved  line)  tends  to  a  logarithmic  for  zUJ  >  0.5  (or  z/h  >  2). 


3.  Results 

3.1.  Time-height  contours 

Figure  4  shows  time-height  colour  contours  of  streamwise  (a)  and  vertical  (b)  velocity 
fluctuations  at  position  x\  (Figure  1(b))  for  case  SCI  (Table  1);  fluctuation  is  defined  as 
deviation  of  a  quantity  from  its  time  average: 

u'(x,t)  =  u(x,t)  —  (u(x,t))t.  (7) 

In  addition,  we  apply  quadrant  analysis  [4,43]  to  the  x  —  z  component  Reynolds  shearing 
stresses  (i.e.  {«',  w'}),  to  illustrate  the  temporal  dynamics  of  turbulent  momentum  fluxes  in 
the  urban  roughness  sublayer.  Following  is  the  established  nomenclature, [43, 44]  we  define 
the  four  quadrants  as 

Q\  (Outward  interaction):  u'  >  0  and  w'  >  0, 

Q2  (Ejection):  u'  <  0  and  u>'  >  0, 

(8) 

Q 3  (Inward  interaction):  u’  <  0  and  u>'  <  0,  and 

Q4  (Sweep):  u’  >  0  and  w'  <  0. 
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Figure4.  Time-height  contours  of  fluctuating  velocity  (u'(x,  t)  =  u(x,  t)  —  (fi(x,  f)),)  components 
for  case  SCI  (Table  1)  at  positions  indicated  in  Figure  1(b):  (a)  u'(xi,  yi,  z,  t)  at  to;  (b)  w'(xi,  yi,  z,  t) 
at  X\,  (c)  fl'(x4,  y4,  z,  t)  at  x4;  and  (d)  u>'(x4,  y4 ,  z,  t )  at  x4.  In  addition  to  colour  floods,  line  contours 
denote  contributions  to  x  —  z  component  of  Reynolds  shearing  stresses,  u'w',  due  to  Q4  ‘sweep’ 
(black)  and  Q2  ‘ejection’  (yellow)  events.  Note  that  in  panels  (a)  and  (b),  u'(x,  y,  z/ H  <  1/8,  f)  =  0 
owing  to  the  presence  of  cubes  of  height  h!H  =  1/8. 


Note  that  in  Figures  4,  5,  and  6,  dimensional  time  has  been  ‘shear  normalised’  by  the 
prescribed  shear  velocity,  uT  and  boundary  layer  depth,  H.  Unless  noted  otherwise,  this 
time  normalisation  is  used  hereafter.  Also,  we  add  that  colour  floods  at  positions  x\  and  x4 
are  shown  since  these  represent  the  limiting  cases  (X2  and  .o,  are  intermediate,  as  seen  in 
Figure  1(b)).  Additional  statistics  below  incorporate  all  positions. 

At  each  vertical  position,  z,  we  determine  the  contribution  from  quadrant  1  to  4  events 
to  the  overall  turbulent  stress  based  on 


{u'w'),'Q(x,,  z ,  t;TL)  =  {u'(xh  z,  t)w'(xh  z,  t)IQ(xh  z;Tl))t, 


(9) 
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Figure5.  Time-height  contours  of  fluctuating  velocity  (u'(x,  t )  =  u(x ,  t)  —  («(x,  f))()  components 
for  case  SC2  (Table  1)  at  positions  indicated  in  Figure  1(b):  (a)  u'(x i,  y\,  z,  t )  atxi;  (b)  w’(x\ ,  vi,  z ,  t) 
at  X\ ;  (c)  u'(x 4,  y4,  z,  f)  at  x4;  and  (d)  w’(X4,  y4,  z,  f)  at  x4.  In  addition  to  colour  floods,  line  contours 
denote  contribution  to  u'w'  due  to  Q4  ‘sweep'  (black)  and  Q2  ‘ejection’  (yellow)  events.  Note  that  in 
panels  (a)  and  (b),  u'(x,  y,  z/H  <  1/4,  t)  =  0  owing  to  the  presence  of  cubes  of  height  h!H=  1/4. 


where  x/  =  local  (two-dimensional)  positions,  x\,  X2,  X3  orx4  (see  Figure  1(b)),  TL  is  the 
so-called  hole  size  and  used  to  specify  threshold,  T,  on  the  magnitude  of  terms  contributing 
to  the  turbulent  stresses,  subscript  Q  denotes  the  quadrant  event  of  interest  (i.e.  sweep  or 
ejection.  Equation  (8)),  and  /g(x7- ,  z;  H)  is  the  indicator  function  which  is  used  to  isolate 
the  role  of  different  events  based  on  their  magnitude  (or  some  predefined  criteria): 


IQ{xj,z-,K)  = 


1  if  |m'(x/,  z)w'(xj,  z)\q  >  T,  and 
0  if  \u'(xj,  z)u>'(xj,  z)\q  <  T, 


(10) 
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Figure  6.  Time-height  contours  of  resolved  and  fluctuating  velocity  for  flow  over  ho¬ 
mogeneous  rough  surface,  case  Ul,  taken  from  centre  of  surface:  (a)  u(xc,  yc,  z,  t )/ur; 
(b)  u'(xc,  yc,  z,  t)/uT  m  u(xc,  yc,  z,  t)/ur  -  {u(xc,  yc,  z,  t))t/uv;  and  (c)  u>'(xc,  yc,  z,  t)/uT  = 
w(xc,  yc,  z ,  t)/ur  —  { w(xc ,  yc,  z,  t))t/uz.  In  addition  to  colour  floods,  line  contours  denote  contri¬ 
bution  to  u'w'  due  to  Q4  ‘sweep’  (black)  and  Q2  ‘ejection’  (yellow)  events. 


where,  in  this  study,  T  =  7Y  =  0,  although  alternative  criteria  can  be  used  to  define  the 
threshold. [44]  By  selecting  T  =  0,  the  contribution  of  all  events  to  ( u'w'),  is  considered. 
This  is  appropriate  for  the  present  purposes,  since  we  wish  to  demonstrate  how  turbulent 
mixing  in  the  roughness  sublayer  varies  with  the  passage  of  coherent,  outer  layer  motions. 
For  boundary  layer  turbulence,  the  ejection  ( Q2 )  and  sweep  (Q4)  events  are  most  relevant 
to  the  shear  stress,  while  the  contributions  from  Q\  and  Q3  events  are  known  to  be  rather 
modest.  [44]  Moreover,  in  the  roughness  sublayer  it  is  known  that  contributions  to  shear 
stresses  associated  with  Q4  events  exceed  Q2  [4,8,10];  in  the  aloft  logarithmic  region,  Q2 
stresses  exceed  Q4.  Thus,  Figure  4(a)  and  4(b)  contains  rich  information  on  the  temporal 
dynamics  of  turbulence  above  cubical  topographies.  First,  Figure  4(a)  («')  shows  a  clear, 
repetitive  advective  lag  in  the  u'  time  series,  illustrated  by  negatively  inclined  low(blue)-  and 
high(red)-momentum  structures.  Thus,  there  is  a  time  lag  between  the  passage  of  coherent 
structures  in  the  outer  layer  and  detection  of  their  effects  in  the  roughness  sublayer,  which 
hereafter  is  called  an  advective  lag.  Figure  4(c)  and  4(d)  show  the  same  flow  quantities 
for  point  X4  (Figure  1(b)).  Flere,  the  11  advective  delay  is  even  more  prevalent,  a  result  of 
coherent  ‘roller’  motions  produced  at  the  cubes  [16]  (.ri)  and  advecting  downwind  [8]  while 
remaining  within  the  sublayer.  We  contend  that  the  onset  of  such  a  momentum  excess  (in 
the  outer  layer)  precedes  a  sublayer  ‘spike’  in  Reynolds  shear  stresses  associated  with  a  Q4 
event.  Rollers  play  an  important  role,  since  they  are  contained  within  the  roughness  sublayer 
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and  facilitate  downward  momentum  fluxes  leading  to  elevated  Reynolds  stresses  owing  to 
sweeps.  Likewise,  the  passage  of  LMRs  in  the  outer  layer  precedes  periods  of  relatively 
subdued  mixing  in  the  roughness  sublayer.  For  example,  consider  the  event  beginning  at 
{tuzH~l  ,z/H}  «  {10.2, 0.78}  and  the  subsequent  advective  lag  before  excitation  at  the  cube 
height,  {tuzH~l  ,z/H}  «  {10.3, 0.1}  (Figure  4(c)).  The  sublayer  excitation  corresponds  with 
significant  Q4  u'w'.  Consultation  of  Figure  4(d)  shows  the  alternating  sign  of  w'  and  the 
congruence  of  this  with  u'u>'.  To  this  extent,  we  note  that  w'  does  not  exhibit  an  advective 
delay.  We  note  also  that  the  sampling  period  duration  for  Figure  4(a)-(b)  is  very  similar, 
tuzH~x  ~  2;  we  showed  different  sampling  periods  for  generality.  This  is  true  also  for 
subsequent  results,  shown  below. 

Finally,  we  emphasise  here  that  periods  of  S’  ^  0  in  the  inertial  layer  (z/H  >  4h/H 
~  0.5)  are  associated  with  the  passage  of  coherent,  quasi-streamwise-elongated  coherent 
motions  [8,25]  in  the  outer  layer.  Thus,  it  is  interesting  to  note  an  apparent  inner-outer 
coupling  between  the  passage  of  such  logarithmic  structures  and  the  excitation  and  subdual 
of  sublayer  mixing.  One  is  reminded  of  the  amplitude  modulation  findings  developed  for 
smooth  wall  turbulent  boundary  layers  regarding  the  ‘imprint’  of  large-scale,  outer  layer 
structures  on  near  wall  motions  [45^18]  (discussion  to  follow  below).  Specifically,  Mathis 
et  al.  [48]  demonstrated  that  inner  (viscous)  layer  streamwise  velocity  fluctuations  are 
subjected  to  an  amplitude  modulation  by  the  passage  of  outer  (logarithmic)  layer  coherent 
structures.  They  developed  a  predictive  model  for  this  amplitude  modulation,  which  required 
inputs  only  from  the  outer  layer  (and  experimentally  determined  empirical  parameters).  The 
results  presented  heretofore  suggest  the  presence  of  an  analogous  amplitude  modulation 
process  for  rough  wall  flows. 

Figure  5  shows  time-height  contours  for  case  SC2  (Table  1 )  at  positions  x\  and  xt, 
(Figure  1(b)).  This  case  is  identical  to  SCI  (Figure  4),  except  that  the  cubes  are  effectively 
doubled  in  size  (i.e.  H/h  =  4  and,  therefore,  L/h  —  8  or  LIH  =  2;  these  case  attributes 
are  summarised  in  Table  1).  Nonetheless,  we  observe  similar  advective  lag  in  the  S'  time- 
height  colour  floods  (Figure  5(b)  and  5(d)).  We  also  observe  that  the  passage  of  a  coherent 
structure,  S'  >  0,  in  the  inertial  layer,  say  z!I /  ss  0.75,  manifests  in  the  roughness  sublayer 
and  canopy  (Figure  5(c))  at  a  later  time.  Likewise,  it  is  evident  that  elevated  Reynolds 
stresses  owing  to  sweeps  at  the  cube  height,  h/H  =  1/4,  are  generally  preceded  by  the  aloft 
passage  of  coherent  u'  >  0  parcels  in  the  logarithmic  layer.  We  observe  similar  w'  patterns 
for  this  case,  with  only  mild  advective  lag  relative  to  the  S' . 

Finally,  it  is  of  interest  to  contrast  Figures  4  and  5  cubical  roughness  results  against 
those  for  case  Ul,  a  homogeneous  rough  wall  (Figure  6).  For  case  Ul,  drag  is  exclusively 
imposed  via  the  equilibrium  logarithmic  law  (Equation  (3)).  We  selected  zo  for  case  Ul  to 
closely  match  the  effective  roughness  length,  z0i  Eff.,  for  case  SC2  (where  z0,  Eft.  is  evaluated 
by  fitting  a  logarithmic  streamwise  velocity  profile  to  the  inertial-layer  velocity  profile). 
Table  2  summarises  zo,Eff.  for  all  cases  in  this  study,  where  zo.Eff.  =  zo  for  cases  Ul  to 
U6.  Thus,  for  case  Ul,  the  roughness  sublayer  and  canopy  processes  responsible  for  mo¬ 
mentum  depletion  are  all  parameterised  by  the  equilibrium  logarithmic  law.  Nevertheless, 
we  report  similar  patterns  in  the  time-height  colour  floods.  Figure  6(a)  shows  the  resolved 
streamwise  velocity,  u,  which  is  included  here  to  demonstrate  that  the  homogeneous  rough 
case  possesses  the  kinds  of  temporal  structural  attributes  observed  for  flow  over  cubes, 
even  though  a  roughness  sublayer  is  not  resolved.  Figure  6(b)  and  6(c)  show  S'  and  w’, 
respectively.  Advective  lag  of  S'  is  apparent  in  the  time-height  contours,  as  is  the  relative 
maxima  in  sweep-driven  Reynolds  stresses  at  the  ‘conclusion’  of  a  S'  >  0  coherent  motion. 
It  is  apparent  for  this  case  also  that  the  u>'  contours  exhibit  negligible  advective  lag. 


Downloaded  by  [Princeton  University]  at  09:53  15  May  2015 


Journal  of  Turbulence 


821 


Figure  7.  Sketch  showing:  (a)  cube  topography  (SCI  to  SC2,  Table  1)  and  (b)  homogeneous 
roughness  transect  and  indication  of  reference  elevation,  zRef.  (U1  to  U6,  Table  1). 
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Figure  8.  Correlation  coefficient,  yuu(xh  z.,  r;zRef )  (Equation  (11)),  for  cases  SC2  (position  x4, 
panel  a)  and  U 1  (panel  b).  Black  markers  denote  shear  normalised  maximum  advective  lag,  rmax  „t  H-i 
(Equation  (12)). 


3.2.  Advective  lag  computation 

We  quantify  the  u'  advective  lag  qualitatively  observed  in  Figures  4-6.  This  is  accomplished 
in  two  steps.  First,  we  must  select  a  reference  height,  zRef ,  at  which  a  reference  data-set  can 
be  collected.  Figure  7  is  a  sketch  of  the  cubic  topography  and  homogeneous  roughness  cases, 
with  indication  of  the  selected  zRef.  positions.  Precise  zRef.  values  are  presented  in  Table  2. 
We  selected  zRef.  to  be  slightly  above  the  ‘top’  of  the  roughness.  We  adopted  this  approach 
since  it  facilitates  comparison  between  these  fundamentally  different  topographies.  From 
here,  we  compute  correlations  maps  as  the  convolutions: 


Yuu(xh  z,  r;zRef.)  =  ( u\xh  z,  t)Iu'(xh  zRef.,  t ))  (r)  = 


z,  t)u'(xi,  ZRef..  t  +  r)dt, 


(ID 
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Figure  9.  Advective  lag  versus  wall-normal  elevation  for  cases  U1  (blue  symbols),  SCI  (black 
symbols),  and  SC2  (red  symbols).  For  cases  SCI  and  SC2,  symbols  correspond  with:  x\  (circles),  X2 
(squares),  x3  (‘plus’  sign),  and  x4  (asterisk).  Blue  symbols  correspond  with  xc  for  the  homogeneous 
rough  case.  Vertical  dashed  lines  denote  ZRef.  for  different  cases  (dashed  blue:  Ul;  dashed  black: 
SCI;  and  dashed  red:  SC2).  Solid  lines  represent  Equation  (17)  predictions  of  time-lag  for  effective 
roughness  lengths  summarised  in  Table  2. 


where  xi  represents  a  spatial  (x  —  y)  position  (Figure  1(b)).  Figure  8(a)  and  8(b)  show 
Yuu(xi,  z,  T;zRef.)  for  cases  SC2  and  Ul,  respectively.  In  Figure  8(a),  the  lower  and  upper 
horizontal  dashed  black  lines  represent  the  cube  height  and  three  times  the  cube  height, 
respectively  (where  three  cube  heights  is  approximately  equal  to  the  roughness  sublayer 
depth).  It  is  clear  that  both  data-sets  exhibit  an  advective  lag,  as  evidenced  by  positive 
correlation  with  aloft  preceding  events.  We  also  observe  similar  patterns  for  the  cubic  and 
homogeneous  rough  cases,  which  is  consistent  with  observations  of  Figures  5(c)  and  6(b). 
We  compute  advective  lag  as 

W.(z;zRef.)  =  argmax((«'(x,,  z,  t)Iu'(xh  zRef.,  0X01  (12) 

t 

where  argmax  denotes  the  value  t  for  which  the  convolution  is  strongest.  We  discard 
rmax.  values  corresponding  with  yuu(xi,  z,  rmax.;zRef.)  <  x>  where  /  =  0.3  is  a  predefined 
threshold.  We  experimented  with  a  range  of  /  values,  finding  generally  that  increasing  the 
threshold  only  serves  to  remove  spurious  values  forz  >  >zRef.  while  the  underlying  rmax(z; 
ZRef  )  trends  were  robust  and  indifferent  to  y .  Figure  9  shows  shear  normalised  advective 
lag,  rmax.(z;  zRei)uTH~l  for  cases  SCI,  SC2,  and  Ul  (case  attributes  found  in  Table  1). 
For  cases  SCI  and  SC2,  it  is  clear  that  positions  x\  to  x4  exhibit  effectively  the  same 
tma profiles.  Above  zRef.,  the  profiles  are  roughly  linear,  and  we  remind  the  reader 
that  ZRef.  is  the  approximate  obstacle  elevation  (or  centre  of  the  mean  streamwise  velocity 
profile  inflection).  Furthermore,  the  advective  lag  is  always  negative  above  zRef ,  which  is 
precisely  consistent  with  qualitative  observations  in  Figures  4-6.  For  position  x\,  there  of 
course  are  no  rmax.Mr/W'  datapoints  below zRef.  due  to  the  solid  cube.  A  few  datapoints  are 
available  for  position  X2,  owing  to  its  position  in  the  lee  of  the  cubes.  Flowever,  at  positions 
X3  and  x4,  a  broad  range  of  rmax  uxH~x  values  are  available.  We  emphasise  also  that  the 
fmax.Mr H~x  values  in  Figure  9  are  only  associated  with  yuu  >  y  (i.e.  they  represent  only 
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Figure  10.  Sketch  of  low-momentum  region  (LMR)  above  array  of  cubes.  The  LMR,  u'  <  0,  is 
denoted  by  transparent  grey.  Encapsulating  the  LMR  are  hairpins  (blue).  The  LMR  is  inclined  at 
angle,  y,  and  the  hairpin  heads  exhibit  positive  spanwise  vorticity  (sketch  denotes  sweep  event,  (34). 
The  individual  hairpin  Tegs’  would  themselves  be  inclined  at  «45°.[24] 


strongly  correlated  datapoints).  Thus,  the  approximate  linearity,  rmax.  ~  — z  forz  >  zRef., 
is  a  product  of  actual  processes  within  the  roughness  sublayer  and  inertial  layer.  In  the 
following  subsections,  we  provide  a  conceptual  explanation  and  scaling  argument  for  the 
underlying  dynamics  responsible  for  the  advective  lag  reported  in  Figure  9. 


3.3.  Passage  of  coherent  motions 

The  presence  of  meandering,  coherent  parcels  of  relatively  low  and  high  momentum  in 
turbulent  wall-bounded  flows  over  smooth  [20-24,26,27]  and  rough  [8,25,50,51]  walls  is 
well  known.  The  LMRs  are  encapsulated  by  hairpin  packets  at  the  interface  between  regions 
of  differing  momentum. [24,28]  For  the  case  of  cube  roughness  such  as  the  cases  considered 
here,  Figure  1 0  is  a  sketch  of  the  aforementioned  structural  attributes.  A  streamwise-vertical 
transect  through  the  LMR  would  reveal  a  typical  inclination  of  y  &  15°. [8, 20, 25, 52]  Quan¬ 
titative  visualisation  of  Figure  10  dynamics  can  also  be  attained  with  vortex  visualisation 
techniques. 

Flere,  we  compute  swirl  strength,  kci,  of  the  instantaneous  three-dimensional  velocity 
field,  u. [20, 53-55]  Swirl  strength  is  computed  by  first  evaluating  the  velocity  gradient 
tensor,  A  —  Vii,  and  computing  its  eigenvalues.  /.„  is  the  imaginary  component  of  the 
second  eigenvalue  of  A.  From  this,  we  obtain  the  signed  swirl  strength,  k*ci,  simply  based 
on  the  vorticity  component  of  interest.  In  the  following,  we  wish  to  study  streamwise- 
wall  normal  vortical  activity,  so  spanwise  vorticity,  coy  =  du/dz  —  dw/dx,  is  the  apropos 
vorticity  component  with  which  to  sign 


A*  =  kc 


COv 


(13) 


Swirl  strength  is  an  illuminating  tool  for  studying  vortical  activity,  since  it  resolves  only  ro¬ 
tation  and  not  shear.  Figure  1 1  shows  contours  of  k*ci  superimposed  on  vector  fields  of  fluctu¬ 
ating  velocity  components,  [w\  w'}.  The  figures  clearly  show  the  presence  of  LMRs  beneath 
a  series  of  hairpin  heads  (evidenced  by  local  regions  of  elevated  A.*).  Figure  1 1(a)  is  for  case 
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Figure  11.  Transects  of  flow  field  statistics  for:  (a)  case  SCI  and  (b)  case  Ul.  Vectors  denote  fluc¬ 
tuating  velocity  components,  {S',  w'},  while  colour  contours  are  signed  swirl  strength,  A*.  (Equation 
(13)).  On  both,  the  solid  black  line  represents  an  approximate  demarcation  between  an  LMR  (below) 
and  an  ambient  fluid  (above)  and  is  inclined  at  y  =  15°.  Local  regions  of  positive  A*  represent  the 
heads  of  hairpins  above  the  LMR,  which  are  followed  by  a  Q2  ‘burst’  event, [49]  as  indicated  above. 


SC2,  and  the  presence  of  immersed  blocks  corresponds  with  regions  of  {u',  w'}  *=»  {0,  0}. 
To  facilitate  discussion,  we  have  added  solid  black  lines  inclined  at  y  =  15°,  which  is  the 
approximate  inclination  angle  of  hairpin  trains  in  such  flows.  [8,50,5 1]  The  coherent  hairpin 
heads  precede  Q2  burst  events  and  we  have  indicated  these  events  for  discussion  on  the 
figures.  Note  also  in  Figure  1 1  (a)  the  presence  of  canopy  vortices  associated  with  separation 
and  cube-scale  coherent  motions,  evidenced  by  elevated  Xci  contours  (for  example,  at  {x/H, 
z/H}  {1.4,  0.15}  and  {3,  0.15}).  For  homogeneous  roughness  case  Ul,  Figure  11(b) 
shows  the  same  quantities  as  Figure  11(a).  Again,  the  presence  of  coherent  LMRs  and 
a  series  of  inclined  hairpin  heads  are  apparent.  We  emphasise  finally  that  the  Figure  1 1 
images  are  taken  from  only  one  timestep  during  LES.  However,  these  images  are  entirely 
representative  of  instantaneous  flow  patterns  observed  at  other  times  during  the  simulation. 

In  Section  3.1  and  Figure  9,  we  showed  vertical  profiles  of  advective  lag,  r 
based  on  a  reference  elevation  just  above  the  canopy  (see  Figures  7  and  9).  The  linearity 
exhibited  by  rmax  wri/_1  points  to  an  underlying  physical  process  in  the  roughness  sublayer 
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and  inertial  layer.  Here,  we  attribute  rmax.  - - z  scaling  to  the  passage  of  regions  of  alter¬ 

nating  low  and  high  momentum  in  the  roughness  sublayer  and  inertial  layer,  as  illustrated 
in  Figure  1 1 . 


3.4.  A  model  for  advective  lag 

Following  the  passage  of  a  LMR,  Figure  1 1  shows  that  the  first  position  in  the  domain  to 
experience  relatively  higher  momentum  would  be  in  the  inertial  layer,  above  the  LMR.  As  the 
LMR  advects  downstream,  relatively  higher  momentum  would  be  recorded  at  progressively 
lower  elevations.  Since  representative  information  is  known  about  the  macroscale  attributes 
of  these  coherent  motions,  we  below  propose  a  simple,  semi-empirical  model  to  predict  the 
advective  lag  between  passage  of  such  motions  in  the  inertial  layer  and  evidence  of  their 
‘imprint’  on  roughness  sublayer  and  canopy  dynamics. 

If  a  LMR  (i.e.  u'  <  0,  above)  has  depth,  S',  its  length  can  be  evaluated  based  on 
assumption  of  a  hairpin  train  inclination  angle: 


Ls  fa  S' /  tan(y) . 


(14) 


Furthermore,  if  we  assume  that  a  representative  advective  velocity  for  LMRs  and  high- 
momentum  regions  (HMRs)  is  the  ‘outer’  velocity,  Uq,  we  can  use  the  equilibrium  loga¬ 
rithmic  law  [40]  and  aerodynamic  roughness  length,  zq,  to  predict  Uq : 


Uo 

U  x 


log 

K 


H 

zo 


(15) 


where  k  is  the  von  Karman  constant.  For  cases  U1  to  U6,  a  large  inertial  region  will 
be  resolved  within  the  computational  domain.  For  cases  SCI  and  SC2,  the  streamwise 
velocity  profile  will  be  inflected  in  the  roughness  sublayer,  before  exhibiting  logarithmic 
wall-normal  scaling  in  the  inertial  layer.  Recall  that  for  cases  SCI  and  SC2,  h!H  =  1/8 
and  1/4,  respectively,  which  may  raise  concern  to  some  readers  about  the  presumption  of 
attaining  a  logarithmic  profile  in  the  computational  domain  (especially  for  SC2).  Figure  3 
shows  mean  streamwise  velocity  profiles  (plane-averaged  and  local)  for  case  SC2,  and  a 
logarithmic  profile  is  clearly  attained  for  z/h  >  2  (Figure  3(b)  shows  the  mean  streamwise 
velocity  profile  with  logarithmic-linear  axis  scaling,  clearly  demonstrating  the  presence  of 
a  logarithmic  profile  in  the  inertial  layer).  We  note  also  Figure  5(a)  from  Coceal  et  al.  [8], 
which  shows  local  streamwise  velocity  profiles  superimposed  on  one  another  and  the  very 
rapid  tendency  to  logarithmic  form  about  the  canopy.  Thus,  an  inertial  layer  exhibiting 
logarithmic  scaling  is  present  for  all  simulations  considered  in  this  work.  rmax  (z;  zRef ) 
is  the  advective  lag  between  passage  of  a  coherent  motion  at  elevation,  z,  and  associated 
modulation  of  processes  at  reference  height,  zRef..  Thus,  Equations  (14)  and  (15)  can  be 
combined  to  obtain  the  advective  lag: 


Anax.  fcZRef.) 


ZRef.  —  Z 
tan(y)  Uq 


^(ZRef.  -  z) 
tan  (y)uT  log (H/zo) 


(16) 
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Figure  12.  Vertical  profiles  of  rmax  (z;  zRef )  for  cases  U2-U6  from  LES  (symbols)  and  Equation 
(17)  predictions  (lines).  Cases:  U2  (blue  squares,  solid  blue  line),  U3  (red  circles,  dashed  red  line), 
U4  (red  squares,  solid  red  line),  U5  (black  circles,  dashed  black  line),  and  U6  (black  squares,  solid 
back  line).  Direction  of  increasing  zq  is  included  for  illustration.  In  the  above,  zRef.///  sa  0.023. 


Normalising  rmax.(z;  zRef.)  by  friction  velocity  and  boundary  depth  (‘shear  normalised’) 
yields 


t|nax.  (z;ZRef.)utH 


-1 


K  (ZRef.  ~  Z) 

tan  (y)H  log(///z0) 


(17) 


Since  we  have  run  the  LES  for  cases  SCI  to  SC2,  Equation  (15)  can  be  used  to  retrieve  an 
appropriate  effective  roughness  length,  zo,  Eft.,  by  matching  the  inertial-layer  profiles  from 
LES  with  predictions  from  Equation  (15).  These  values  are  summarised  in  Table  22 .  Follow¬ 
ing  deduction  of  zo.Etr.-  Uq  is  recorded  also  for  substitution  in  Equation  (17).  Finally,  we  use 
Y  =  15  deg  for  Equation  (17)  predictions  considered  here,  which  is  comparable  to  values 
for  rough  wall  flows. [8, 25, 50-52]  For  cases  SCI,  SC2,  and  Ul,  Figure  9  shows  predictions 
from  Equation  (17)  (solid  lines).  For  zRef JH  <  z/H  <  0.2,  Equation  (17)  predictions  exhibit 
reasonable  agreement  with  the  LES  results  for  case  U 1 .  This  is  the  case  of  homogeneous 
roughness  in  which  the  roughness  sublayer  may  be  z/H  <  0.1.  For  cases  SCI  and  SC2, 
Equation  ( 17)  predictions  agree  somewhat  well  with  the  LES  data  for  zRef.///  <  z/H  <  0.4, 
corresponding  with  2  to  4  times  the  cube  height.  Thus,  the  model  performance  is  moderately 
successful  in  the  roughness  sublayer  and  into  the  inertial  layer.  To  broaden  the  test  cases, 
we  include  Figure  12,  which  shows  rmax.(z;  zRcr.)  from  cases  U2  to  U6  and  predictions 
from  Equation  (17).  Once  more,  we  observe  moderate  agreement  between  the  model  and 
LES  data  for  zRef JH  <  z/H  <  0.2.  Note  also  the  apparent  monotonic  increase  in  slope  of 
the  advective  lag  profile,  3rmax.(z;  zRef.)/3z,  with  increasing  aerodynamic  roughness  length 
(partial  derivative  here  is  a  reference  to  assertions  that  aspects  of  the  Equations  (14)— (17) 
development  could  be  generalised,  for  example  the  choice  of  advective  velocity  and  incli¬ 
nation  angle  [22]).  We  attribute  this  to  intensification  of  roughness  sublayer  processes  with 
greater  imposed  drag.  To  this  extent,  we  note  that  3rmax.(z;  zRef.)/3z  for  sase  SC2  exceeds 
sase  SCI  (see  Figure  9);  the  cube  height  for  sase  SC2  is  double  that  of  SCI.  Note  also 
that  zRef.///  ~  0.023  for  sases  U2-U6  (Figure  12),  although  the  nominated  z0j  E sJH  varies 


Downloaded  by  [Princeton  University]  at  09:53  15  May  2015 


Journal  of  Turbulence 


827 


Figure  13.  Vertical  profiles  of  rm x.(z;  zRef.)  for  flow  over  topography  SC2,  resolved  with  Nx  = 
Nv  =  N:  =  64  (light  grey  symbols),  Nx  =  Nv  =  Nz  =  128  (dark  grey  symbols),  and  Nx  =  Nv  = 
Nz  =  256  (black  symbols).  Solid  black  line  is  Equation  (17)  prediction  with  z$IH  =1.08  x  10~2 
(Table  2).  Symbols  correspond  with  position  .to  ( +  ),  position  x2  (❖),  position  x3  (o),  position  x4  (*) 
in  Figure  1(b). 


considerably.  Increasing  roj  mJH  is  tantamount  to  ‘more  rough’  topographies  with  greater 
heights. 


3.5.  Resolution  testing 

A  final  component  of  this  study  involves  evaluating  the  extent  to  which  resolution  -  spatial 
and  temporal  -  influences  the  advective  lag  profiles.  In  the  results  presented  so  far,  case 
SC2  (Table  1)  was  resolved  with  Nx  =  Ny  =  Nz  =  128.  For  generality,  we  have  also 
modelled  flow  over  the  SC2  topography  with  Nx  =  Ny  =  Nz  —  64  and  Nx  =  Nv  =  Nz  = 
256.  The  resulting  advective  lag  profiles  are  shown  in  Figure  13,  where  the  symbol  colour 
corresponds  with  different  resolutions.  Included  on  this  figure  is  the  Equation  (17)  profile 
for  the  appropriate  effective  roughness  length.  It  is  apparent  that  the  advective  lag  profiles 
(symbols)  and  model  profile  (solid  line)  exhibit  the  same  slope.  Moreover,  no  significant 
resolution  dependence  is  apparent.  As  per  the  procedure  outlined  in  Section  3.2,  we  omit 
advective  lag  computations  for  cases  in  which  yuu  <  /. 


4.  Conclusion 

LES  with  an  IBM  has  been  used  to  model  neutrally  stratified  atmospheric  boundary  layer 
flows  over  urban-like  topographies.  We  also  consider  a  suite  of  cases  with  differing  ho¬ 
mogeneous  momentum  roughness,  without  the  presence  of  urban-like  cubes.  We  recorded 
time-series  fluctuating  velocity  statistics  across  the  depth  of  the  boundary  layer,  providing 
data-sets  that  would  typically  emerge  from  field  studies  using  a  micrometeorological  tower 
and  series  of  sonic  anemometers.  Flowever,  owing  to  the  use  of  LES,  we  can  achieve  high 
resolution  measurements  in  both  time  and  elevation  (at  predefined  locations).  With  this, 
we  considered  contour  maps  of  fluctuating  streamwise  and  vertical  velocity  components 
with  respect  to  shear-normalised  time  and  elevation;  superimposed  on  these  images  were 
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contours  of  Reynolds  shearing  stresses  owing  to  turbulent  ‘sweeps’  and  ‘ejections ’.[43] 
The  colour  flood  contours  clearly  revealed  the  presence  of  an  advective  lag  between  mo¬ 
mentum  excess!  deficit)  in  the  inertial  layer  and  excitation(subdual)  of  flow  processes  in  the 
roughness  sublayer.  This  advective  lag  was  observed  for  all  cases  -  cube  topography  and 
homogeneous  roughness. 

We  computed  the  advective  lag,  rmaxuzH~l,  and  observed  a  linear  scaling  against 
elevation.  This  prompted  consideration  of  underlying  physics  responsible  for  the  advective 
lag.  We  used  spatial  flow  and  vortex  visualisations  to  argue  that  coherent  outer  layer  motions 
(LMRs  and  HMRs)  are  responsible  for  the  advective  lag  between  sublayer  and  inertial-layer 
dynamics.  With  this,  we  used  canonical  attributes  of  the  LMRs  (hairpin  head  inclination 
angle,  y,  and  advective  velocity,  Uq),  to  develop  scaling  arguments  for  the  advective  lag. 
We  observed  promising  results. 

This  work  is  inspired  by  the  amplitude  modulation  concepts  developed  in  recent  years 
by  Marusic  and  co-workers  [45-48],  The  existing  literature  on  amplitude  modulation  is 
in  the  context  of  smooth  wall  flows  (thus,  ‘inner’  refers  to  the  viscous  sublayer,  not  the 
roughness  sublayer).  Nonetheless,  it  is  interesting  to  observe  some  qualitative  similarities  in 
the  present  findings.  Moreover,  the  present  work  has  been  performed  in  the  context  of  urban- 
like  topographies,  although  one  could  imagine  applications  to,  for  example,  control  of  very 
large  wind  farms  over  which  the  atmospheric  boundary  layer  is  fully  developed.  Indeed, 
the  modelling  framework  (Equation  (17))  presented  herein  requires  only  two  inputs:  y  and 
zo.  Ref.-  We  used  a  somewhat  standard  value  of  y,  and  computed z0j  Ref.  by  fitting  logarithmic 
velocity  profiles  to  the  time-  and  plane-averaged  streamwise  velocity  profiles.  We  expect  that 
for  urban  topographies  exhibiting  greater  height  variability  than  the  present  (uniform  height) 
cases,  or  for  entirely  different  classes  of  complex  topography,  the  modelling  arguments 
presented  here  would  be  valid. 
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Notes 

1.  Throughout  this  study,  we  adopt  the  following  nomenclature:  u  =  [ii,  v,  wj  is  the  streamwise, 
spanwise,  and  vertical  velocity  component,  respectively;  x  =  [x,  y,  z)  is  streamwise,  spanwise, 
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and  vertical  position,  respectively;  and  indices  (where  used)  i  =  1,2,  and  3  are  streamwise, 
spanwise,  and  vertical  direction,  respectively. 

2.  Alternatively,  one  could  use  predictive  models  for  z0  based  on  attributes  of  the  topography 
[56-58]  although  this  would  require  empirical  parameters. [59] 
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